###############################
# Analysis of conjoint data 
###############################

rm(list=ls())

# load packages 
library(foreign)
library(lmtest)
library(sandwich)
library(stargazer)
library(msm)

sink("~/Dropbox/Atacama 2015/10_replication/004_conjoint_diagnostic.txt")

# load function that does clustered SEs
vcovCluster <- function(
  model,
  cluster
)
{
  require(sandwich)
  require(lmtest)
  if(nrow(model.matrix(model))!=length(cluster)){
    stop("check your data: cluster variable has different N than model")
  }
  M <- length(unique(cluster))
  N <- length(cluster)           
  K <- model$rank   
  if(M<50){
    warning("Fewer than 50 clusters, variances may be unreliable (could try block bootstrap instead).")
  }
  dfc <- (M/(M - 1)) * ((N - 1)/(N - K))
  uj  <- apply(estfun(model), 2, function(x) tapply(x, cluster, sum));
  rcse.cov <- dfc * sandwich(model, meat = crossprod(uj)/N)
  return(rcse.cov)
}

##############################
# Prepare data
##############################

# load data
d = read.dta("~/Dropbox/Atacama 2015/10_replication/conjoint_paipote_final.dta")

dim(d)

# drop failed conjoints
d = d[d$failcon!=1,]

# gen treatment indicator
d$t_ind = NA
d$t_ind[d$a2>2]=1
d$t_ind[d$a2<3]=0
table(d$t_ind)

# gen control indicator
d$c_ind = 1 - d$t_ind 
table(d$c_ind)

# gen affected area indicator 
table(d$zone)
d$zone2 = NA
d$zone2[d$zone==1] = 1
d$zone2[d$zone==3] = 1
d$zone2[d$zone==2] = 0
table(d$zone2)

# pairs
d_pair1 = d[d$pair==1,]
d_pair2 = d[d$pair==2,]
d_pair3 = d[d$pair==3,]
d_pair4 = d[d$pair==4,]
d_pair5 = d[d$pair==5,]
d_pair6 = d[d$pair==6,]

######################################################
## ACME using regression, center as reference category
######################################################

# Check levels
print(levels(d$atideology))  
d$atideology = factor(d$atideology,levels(d$atideology)[c(2,1,4,3)])
print(levels(d$atideology)) 

##################
# Diagnostics
##################

# Carryover effects
d$fpair = as.factor(d$pair)
d1 <- lm(selected ~ atideology + atprofession + atgender + atage + atexperience + atexpectations + 
                    fpair + 
                    atideology*fpair + atprofession*fpair + atgender*fpair + atage*fpair + atexperience*fpair + atexpectations*fpair,
                    data=d)
d1_c <- round(coeftest(d1, vcov = vcovCluster(d1, cluster = d$idnum)),4)
d1_c 

# Random assigment: regress gender on attributes (regress covariate on treatment)

# gender
d2 <- lm(c1 ~ as.factor(atideology) + as.factor(atprofession) + as.factor(atgender) + as.factor(atage) + as.factor(atexperience) + as.factor(atexpectations), data=d)
d2_c <- round(coeftest(d2, vcov = vcovCluster(d2, cluster = d$idnum)),4)
d2_c 

# Profile order effects
d$fcandidate = as.factor(d$candidate)
d5 <- lm(selected ~ atideology + atprofession + atgender + atage + atexperience + atexpectations + 
                    fcandidate + 
                    atideology*fcandidate + atprofession*fcandidate + atgender*fcandidate + atage*fcandidate + atexperience*fcandidate + atexpectations*fcandidate,
                    data=d)
d5_c <- round(coeftest(d5, vcov = vcovCluster(d5, cluster = d$idnum)),4)
d5_c

######################
# Regression tables
######################

# Random assigment:gender 
d2_c
stargazer(d2_c,title="Balance test: Gender",
          type = "text",
          align=TRUE, 
          omit.stat=c("LL","ser","f"), 
          covariate.labels=c("Right","Left","Independent","Teacher","Engineer","Female","40","50","Council member","Mayor","Will distribute a financial relief"),
          no.space=TRUE,  
          table.placement = "H",
          #single.row = TRUE,
          star.cutoffs = c(0.05,0.01,0.001),
          dep.var.caption  = "Outcome",
          dep.var.labels   = "Gender")
          
# Candidate order effects
d5_c
stargazer(d5_c,title="Candidate order effects",
          type = "text",
          align=TRUE, 
          omit.stat=c("LL","ser","f"), 
          keep = c(13:23),         
          covariate.labels=c("Right*Candidate 2",
                             "Left*Candidate 2",
                             "Independent*Candidate 2",
                             "Teacher*Candidate 2",
                             "Engineer*Candidate 2",
                             "Female*Candidate 2",
                             "40*Candidate 2",
                             "50*Candidate 2",
                             "Council member*Candidate 2",
                             "Mayor*Candidate 2",
                             "Will distribute a financial relief*Candidate 2"), 
          no.space=TRUE,  
          table.placement = "H",
          #single.row = TRUE,
          star.cutoffs = c(0.05,0.01,0.001),
          dep.var.caption  = "Outcome",
          dep.var.labels   = "Electoral Choice")     
          
# Carryover effects          
d1_c    
stargazer(d1_c,title="Profile order effects",
          type = "text",
          align=TRUE, 
          omit.stat=c("LL","ser","f"), 
          keep = c(19:39),     
          covariate.labels=c("Right*Pair 2","Left*Pair 2","Independent*Pair 2",
                             "Right*Pair 3","Left*Pair 3","Independent*Pair 3",
                             "Right*Pair 4","Left*Pair 4","Independent*Pair 4",
                             "Right*Pair 5","Left*Pair 5","Independent*Pair 5",
                             "Right*Pair 6","Left*Pair 6","Independent*Pair 6",
                             "Right*Pair 7","Left*Pair 7","Independent*Pair 7",
                             "Right*Pair 8","Left*Pair 8","Independent*Pair 8"),     
          no.space=TRUE,  
          table.placement = "H",
          #single.row = TRUE,
          star.cutoffs = c(0.05,0.01,0.001),
          dep.var.caption  = "Outcome",
          dep.var.labels   = "Electoral Choice")  

sink()
          